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We characterize the distributions of size and duration of avalanches propagating in complex networks. By 
an avalanche we mean the sequence of events initiated by the externally stimulated 'excitation' of a network 
node, which may, with some probability, then stimulate subsequent firings of the nodes to which it is connected, 
resulting in a cascade of firings. This type of process is relevant to a wide variety of situations, including neuro- 
science, cascading failures on electrical power grids, and epidemology. We find that the statistics of avalanches 
can be characterized in terms of the largest eigenvalue and corresponding eigenvector of an appropriate adja- 
cency matrix which encodes the structure of the network. By using mean-field analyses, previous studies of 
avalanches in networks have not considered the effect of network structure on the distribution of size and du- 
ration of avalanches. Our results apply to individual networks (rather than network ensembles) and provide 
expressions for the distributions of size and duration of avalanches starting at particular nodes in the network. 
These findings might find application in the analysis of branching processes in networks, such as cascading 
power grid failures and critical brain dynamics. In particular, our results show that some experimental signa- 
tures of critical brain dynamics (i.e., power-law distributions of size and duration of neuronal avalanches), are 
robust to complex underlying network topologies. 

PACS numbers: 



I. INTRODUCTION 

In this paper we study the statistics of avalanches propagat- 
ing in complex networks. The study of avalanches of activity 
in complex networks is relevant to diverse fields, including 
epidemiology 01 Q], genealogy yD, and neuroscience 0-fHtl- 
The sim ples t case of an avalanche corresponds to a branching 
process Il4ll5ll . first studied by Galton and Watson [3], which 
can be considered as an avalanche propagating in a tree net- 
work. Various generalizations to the case where avalanches 
propagate in a more general network have been considered 
recently lfl3l[T6|jl8ll . and related problems such as the distri- 
bution of cluster size in percolation models lfl 9[[2 0ll and self- 
organized criticality in the "sandpile" model |21] have been 
studied. In contrast to these previous studies, we develop a 
theory of avalanche size and duration on complex networks 
that, instead of using some form of mean field analysis, ex- 
plicitly includes the network topology. This approach allows 
for an analysis of avalanches starting from arbitrary nodes in 
the network and the effect of nontrivial network structure on 
the distribution of avalanche size and duration. 

Our formalism in this paper is general, describing dynam- 
ics with applications to a wide variety of systems. Our results 
are correspondingly general, but they may be of particular in- 
terest to those investigating recent experimental observations 
of avalanches of neuronal bursting in the mammalian cortex. 
When a neuron fires, it stimulates other neurons which may 
subsequently fire. When this linked activity occurs in a cas- 
cade, it is called a neuronal avalanche (experimentally, neu- 
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ronal avalanches are observed propagating in functional net- 
works where each node represents a group of neurons). Re- 
cent experiments have studied neuronal avalanches of activity 
in the brains of awake monkeys fl, anesthetized rats rf22ll . 
slices of rat cortex [HI, and humans 1E3I . These stud- 
ies found that when the tissue is allowed to grow and oper- 
ate undisturbed in homeostasis [6], both the size and temporal 
duration of neuronal avalanches is distributed according to a 
power-law. In contrast, the application of drugs that selec- 
tively decrease the activity of inhibitory or excitatory neurons 
results in avalanches with different statistics ]5|]. Based on 
these observations, it has been argued and demonstrated ex- 
perimentally that many neuronal networks operate in a crit- 
ical regime that leads to power-law avalanche distributions 
B [1111221 l23ll . maximized dynamic range JH 0-191], an ^ max- 
imized information capacity lflol - [l2ll . Therefore, it is of great 
interest to characterize this critical state and to understand 
how experimental signatures of criticality may change upon 
modification of the underlying network (e.g., changes induced 
by the drugs used in experiments). 

We find that the statistical properties of avalanches are de- 
termined by spectral properties of the matrix whose entries 
A mn are the probabilities that the avalanche propagates from 
node n to node m. In particular, the eigenvalue A of max- 
imum magnitude (by the Perron-Frobenius theorem A is real 
and positive if A mn > 0) and its associated eigenvector play a 
prominent role in determining the functional form and the pa- 
rameters for the statistical distribution of avalanche size and 
duration. While many of our findings have analogous results 
in classical Galton- Watson branching processes lfl4l Il5ll . we 
emphasize that our analysis allows us to identify how changes 
in network structure affect the parameters of the statistical dis- 
tributions of avalanche size and duration. Moreover, our the- 
ory allows us to obtain the statistics of avalanches starting at 



particular network nodes. 

This paper is organized as follows. In Sec. |II]we describe 
our model for avalanche propagation in networks. In Sees. [HI] 
and [IV] we analyze the statistics of avalanche duration and 
size. In Sec. [V] we validate our analysis through numerical 
experiments. Section [VT]presents further discussion and con- 
clusions. 



II. FORMULATION 



To model the propagation of avalanches in a network, we 
consider a network of N nodes labeled m = 1, 2, N, Each 
node m has a state x m = or 1. We refer to x m — as the 
resting state and to x m — 1 as the excited state. At discrete 
times t = 0, 1, the states of the nodes x^ are simultane- 
ously updated as follows: (i) If node m is in the resting state, 
i* n = 0, it can be excited by an excited node n, x^ = 1, with 
probability < A mn < 1, so that x 1 ^ 1 = 1. (ii) The nodes 
that are excited, x* n = 1, will deterministically return to the 
resting state in the next time step, = 0. We therefore 

describe a network of N nodes with a N x N weighted net- 
work adjacency matrix A = {A mn }, where A mn > may 
be thought of as the strength of connection from node n to 
node m, and A mn = implies that node n does not connect 
to node m. We will assume that given any two nodes n and 
m, the probability that an excitation originating at node n is 
able to excite node m (through potentially many intermediate 
nodes) is not zero. This is equivalent to saying the network is 
fully connected, and therefore the matrix A is irreducible. 

Starting from a single excited node k (x° = 1 if m = n 
and x® n = if m ^ n), we let the system evolve according to 
the dynamics above, and observe the cascade of activity until 
there are no more excited nodes. This motivates the following 
definitions, which are illustrated in Fig. Q]: (1) an avalanche is 
the sequence of excitations produced by a single excited node; 
(2) the duration d of an avalanche is defined as the total num- 
ber of time steps spanned by the avalanche: if the avalanche 
starts with x® = 1, then 



minis; 



k 



Oforallfc}. 



(1) 



An avalanche that continues indefinitely is said to have infi- 
nite duration; (3) the size x of an avalanche starting at x^ = 1 
is defined as the total number of nodes excited during an 
avalanche, allowing for nodes to be excited multiple times: 



d-l N 

EE 4. 

t=0 k=l 



(2) 



Note that it is possible for an avalanche to have size larger 
than the total size of the network (e.g., if d n = oo, then x n = 
oo). Our goal in this paper is to determine the probability 
distributions of these variables in terms of the matrix A. 




FIG. 1: An example avalanche is shown, where circles represent 
nodes, arrows represent links, and numbers inside nodes correspond 
to the time step at which each node is activated. Starting from a sin- 
gle excited node, labeled 1, the avalanche spreads to two other nodes, 
labeled 2, and so on. Note that the presence of a link does not guar- 
antee the transmission of excitation. The example avalanche above 
lasts for five time steps and excited a total of six nodes in addition to 
the initial node, so d — 5 and x = 7. 



III. DISTRIBUTION OF AVALANCHE DURATION 

In order to analyze the statistics of avalanche duration, we 
define c n (t) as the probability that an avalanche starting at 
node n has duration less than or equal to t, 



C(t) = P(dn < *)• 



(3) 



The quantity c n (t) is the cumulative distribution function 
(CDF) of the random variable d n . In what follows, we will 
restrict our attention to a class of networks that we call lo- 
cally tree-like. By locally tree-like, in this paper we shall 
mean that, for any given t not too large, and pair of nodes 
j and k, if there exists a directed path of length t from j to 
k, then it is rarely the case that there will also exist a second 
such path l24ll . Many networks found in applications are of 
this type, and it has been found that the locally tree-like ap- 
proximation works very well in describing various dynamical 
processes while still capturing the effects of network hetero- 
geneity JH I24U26I1 . For these networks, we can approxi- 
mately treat the avalanches propagating to different neighbors 
of node n as independent, and write the recursion relation 



c„(t+l) = 



n[< 

m— 1 



1 — A mn ) + A r 



(4) 



together with c„(0) = which follows from the definition (0. 
The right hand side of Eq. (|4]i is the probability that nodes are 
either not excited by node n, or, if they are, that they generate 
avalanches of duration at most t: (1 — A mn ) is the probabil- 
ity that an excitation does not pass from node n to node m, 
whereas A mn c m (t) is the probability that an excitation does 
pass from node n to node m and the resulting avalanche has 
duration at most t. Note that Eq. (01 can treat any node n 
as the starting node for an avalanche. As discussed above, 
Eq. assumes that the descendent branches of the avalanche 
are independent. It is, however, possible that an avalanche 



3 



may branch in such a way that two branches interact at a later 
time. Nevertheless, for the networks we studied we found that, 
while these events do occur for large avalanches, they do not 
significantly affect our predictions. We show numerical re- 
sults confirming this in Sec. [V] 

We are interested in the distribution of long avalanche du- 
ration, i.e., in the asymptotic form of c n (t) for t — > oo. By 
definition (see also AppendixlAt. c n (t) is a bounded, increas- 
ing function of t, and therefore it must converge to a value 
lim^oo c n (t) = b n < 1 which can be interpreted as the 
probability that an avalanche starting at node n has finite dura- 
tion. Our analysis will be based on whether or not this limit is 
strictly less than one or equal to one. As shown in Appendix 
lAl this is determined by the Perron-Frobenius eigenvalue of 
A, A: if A < 1, then lim^oo c n (t) = 1. The case A < 1 will 
be referred to as the subcritical case, and the case A = 1 will 
be referred to as the critical case. On the other hand, if A > 1, 
then lim^oo c n (t) = b n < 1, which implies that there is 
a nonzero probability that an avalanche has infinite duration. 
This case will be referred to as the supercritical case. The 
asymptotic form of c n (t) will be analyzed separately for these 
three cases below. 



A. Subcritical Networks (A < 1) 

In the subcritical case, b n = 1 is the only fixed point of the 
system Eq. © (see Appendix [A}. To analyze the asymptotic 
form of c n (t), we assume it is close to the fixed point and de- 
fine the small quantity f„ (t) = 1 — c n (t). Linearizing Eq. © 
we obtain 



N 



(5) 



m — 1 



Assuming exponential decay (or growth) of perturbations, 

fn{t) = X t v n , we obtain 



N 



Xv n = 2J A r, 



(6) 



1 is linearly stable when A < 1. 

The probability density function (PDF) of avalanche dura- 
tion is given by p„(t) = P(d n = t) = c n (t) - c n (t - 1), 
so 



P. 



,{t) ~ (A- 1 - l)v n X\ 



(8) 



which decays exponentially to zero with decay rate ln(l/A). 

In summary, we can draw two predictions from the analy- 
sis above for subcritical networks: (i) the PDF of avalanche 
duration decays exponentially towards zero as A*, and (ii) the 
probability that an avalanche started at node n lasts t steps is 
proportional to the n th entry of the left eigenvector of A, v n . 
These predictions are tested in Sec. [V] 



B. Supercritical networks (A > 1) 

A linear stability analysis of the fixed point b n — 1 in the 
supercritical case shows that this fixed point is linearly un- 
stable. This implies (see Appendix lAl that there exists an- 
other fixed point b n to which c n (t) converges from below, 
lim t _ i . 00 c„ (t) = b n < 1. Thus, there is a nonzero probability 
that an avalanche will have infinite duration. Our analysis be- 
low characterizes the distribution of finite avalanche duration 
in supercritical networks. We first note that the fixed point b n 
satisfies 



N 



|^(1 A mn ) + A mn b ri 



(9) 



Again, we introduce the quantity f n (t) = b n — c n (t), and con- 
sider the limit when t is large and /„ is small. We substitute 
this into Eq. © and rewrite it as 



■/n(* + l) 



TV 



(1 A mn j -\- Auinb^ 
By defining a new matrix D with entries 



(10) 



Thus, A is an eigenvalue of A and v = [v±, V2, vjv] its left 
eigenvector. We identify A as the Perron-Frobenius eigenvalue 
since, having the largest magnitude among all the eigenvalues, 
X t v n will be the dominant term as t — > oo when compared 
with the other modes. We note that for finite t, this approx- 
imation is good as long as there is a large enough separation 
between A and the rest of the spectrum of A. This issue is 
discussed in ll27ll . where it is found that this separation is typ- 
ically large in networks without strong community structure. 
Henceforth, we will assume that A is well separated from the 
rest of the spectrum of A. Therefore, c„ (t) approaches 1 ex- 
ponentially as 



Cn{t) W 1 



X t V r , 



(7) 



where v is the left eigenvector of A corresponding to A; v n > 
by the Perron-Frobenius theorem |J28| . The fixed point b n = 



(1 A mn ^ -\- A ran b Tl 
and linearizing Eq. (TTOb we find 



(ID 



(12) 



As in the subcritical case, we conclude that f n (t) re A^,w„, 
where w is the left Perron-Frobenius eigenvector of the ma- 
trix D and Ad its Perron-Frobenius eigenvalue. As argued in 
AppendixlBl A^ < 1 when A > 1, thus ensuring exponential 
convergence. Therefore, we have 

Cn{t) w b n -WnXp. (13) 

As in the subcritical case, the probability density function 
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(PDF) of avalanche duration is given by 

Pn(t) ~ (A^ 1 - l)u>„A^, 



(14) 



which decays exponentially to zero with decay rate 1ii(1/Ad). 

In summary, for supercritical networks: (i) the PDF of 
avalanche duration decays exponentially towards zero as A^,, 
and (ii) the probability that an avalanche started at node n lasts 
t steps is proportional to the n th entry of the left eigenvector 
of D, w n - These predictions are tested in Sec. [V] We note 
that these predictions simplify to those drawn from Eq. (|8) 
if the network is subcritical, in which case b n = 1, Eq. ( fTTT i 
simplifies to D mn = A mn , and therefore \ D = A and w = v. 



C. Critical Networks (A = 1) 

The analyses above show that if A = 1, the fixed point 
b n = 1 is marginally stable. This fixed point must be an at- 
tracting fixed point, since c n (t) is nondecreasing and b n = 1 
is the only fixed point of Eq. (HJi as shown in Appendix [A] 
To determine the asymptotic form of c n (t) for large t, we 
let c n (t) = 1 — f n (t). We assume that Eq. (01 has a solu- 
tion whose asymptotic functional form in t (to be determined) 
can be extended to a differentiable function of a continuous 
time variable t. Since the convergence of /„ (t) to is slower 
than exponential, we look for a solution /„(<) which is slowly 
varying in t when /„ (t) is small, and approximate 



/n(*+l)«/n(t) + 



(15) 



The slowly varying assumption implies that df n (t)/dt = 
fn(t) *C f n (t) as /„ — > 0, which we note excludes an ex- 
ponential solution. Substituting Eq. (fT~5T > into Eq. ©, we get 



N 



(16) 



Assuming f n (t) <C 1 and expanding to second order, we get 
after simplifying and dropping the time notation for clarity, 

fn + f n ~ ^ ] A mn f m — — y ' y ] AmnAknfmfk- (17) 

m m k^m 

The leading order terms are /„ on the left-hand side and 
A mn f m on the right-hand side, so for these to balance 
as / — > requires 



fn — y ] A n 



,nfn 



(18) 



error from this limit solution, so we set 

f n {t) = K{t)v n /{v)+e n {t), 



(19) 



where we assume e n fn(t)> £ ' n ^ /n(*)> an< ^ me 
term (v) = J2n=i V n/N is included to make K{i) in- 
dependent of the normalization of v. Inserting this in 
Eq. (fTTI i, neglecting terms of order e 1 , e 2 , fe, and using 

E m Lfe^m A mn A kn v m v k k vl, we obtain 



JV 



e n + K'{t)v n /(v) = A mne m - ^K 2 (t)v 2 J{v) 2 . (20) 



To eliminate the unknown error term e, we multiply by u n , 
where u is the right eigenvector of A satisfying Au = u, and 
sum over n. The error terms cancel and we obtain an ordinary 
differential equation (ODE), 



K'(t) = 



1 (uv 2 ) 

2 (uv){v) 



K 2 (t), 



where (xy) = ^ J2 n x nUn- Solving this ODE yields 



K(t) 



P ^ 2 (uv){v) L 



(21) 



(22) 



where f3 is an integration constant. In terms of the original 
variables, we obtain 



C(t) w 1 



+ iJS^Krt' 

r 2 {uv} {v} 



(23) 



The PDF, in the continuous time approximation, is given by 

Pn(t) = C' n (t), 



p n (t) oc 



P ^ 2 (uv)(v) L 



(24) 



From Eq. (l24l we make the prediction that as t — > oo, 
Pn{t) ~ v n t~ 2 - This prediction is tested in Sec. [V] 



IV. DISTRIBUTION OF AVALANCHE SIZE 

In order to analyze the distribution of avalanche size, we de- 
fine the random variable x n as the size of an avalanche starting 
at node n. Let z mn be a random variable which is 1 if node n 
excites node m and otherwise, so that z mn = 1 with proba- 
bility A mn and with probability 1 — A mn . Thus 



Therefore, in this limit the vector f(t) = 
[fi(t), /2 (<),..., Jn {t)] T has to be proportional to the 
normalized left eigenvector v of A with eigenvalue A = 1. 
Thus, a slowly varying solution only exists for a critical 
network. Since v is independent of time, the constant of 
proportionality must be time dependent, f n (t) = K{t)v n . 
Now, for finite /, we expect the solution to deviate by a small 



N 



Xn = 1 + X! 



(25) 



When A > 1 there is a nonzero probability that an avalanche 
has infinite duration, and therefore infinite size, as demon- 
strated in Sec. IIII Bl and Appendix lAl Therefore, we will re- 
strict our attention only to the distribution of avalanches that 
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are finite. To study this distribution, we define the moment 
generating function 



(f> n (s) = E[e~ 



\X n < OO 



(26) 



We now use Eq. (f25b to derive a relation between the mo- 
ment generating functions corresponding to different nodes. 
First, we rewrite the condition x n < oo for node n in terms 
of events applicable to its neighbors. An avalanche starting 
at node n is finite if and only if for every node to, either (i) 
the excitation does not pass from node n to node to, or (ii) 
the excitation passes from node n to node to but the subse- 
quent avalanche starting from node m is finite. Therefore, 
we rewrite the condition x n < oo as the requirement that 
for any to, (z mn ,x m ) G Z mn U W mn , where we have de- 
fined the disjoint sets of events Z mn = {z mn = 0} and 
Wmn = {im < oo and z mn = !}■ Assuming the inde- 
pendence of the random variables x m (consistent with the lo- 
cally tree-like assumption used in the previous section), we 
can rewrite <fi n (s) as 



N 

n4 



U W n 



rn — 1 



(27) 



where the expectation E[] is taken over realizations of the 
random pairs (z mn ,x m ). Denoting P(W) as the probability 
of an event set W, we relate the expected value in the product 
in Eq. d27b to the probabilities of the events W mn and Z mn : 

E[e~ SZm " Xm \Z mn U W mn }P(Z mn U W mn ) 
= E[e- sz ™ x ™\Zmn]P{Zmn) 
+ E[e- sz ™ x ™ \Wmn]P{W mn ). (28) 

Using the following relations that follow from the definitions 
above, 

P '{Z mn ) 1 Amni 
P\Wmn) — A mn b mi 
P(Z mn U W mn ) = (1 - Amn) + A 
E[e- sz ™ x ™\Wmn] = 4>m{s), 

B[e _ " Mt "W = l 1 



Substitution into Eq. ( |28i > gives 



E[e~ sz ™ Xm \Z mn UW mn ] 

(1 A mn j -\- b rn A rnri (j) rn 

(s) 



(29) 



(1 A mn ) -\- h Tri A Trin 
Inserting this into Eq. (|27T i we obtain one of our main results, 



N 



bn(s)=e-» J] K — 



(l A mn j -\- 6 m j4 mn ^) m (s) 



(30) 



(1 Amn ) 4~ bmAjmi 

Defining g n (s) — <fi n (s) — 1, and the matrix H with entries 



bm.A T 



(1 A rnn ) + b m A Ti 
we can rewrite Eq. fl30b as 



N 



1 + 9n(s) 



m— 1 



+ i?mnffm(s)]. 



(3D 



(32) 



Defining the N xN matrix, B = diag(6i, 62, 6jv), we have 
from Eqs. O and (EB, that HE,- 1 = B^D. Thus the ma- 
trix H is related to the matrix D by a similarity transformation 
and thus has the same spectrum. Therefore, we will denote the 
Perron-Frobenius eigenvalue of H by Ad . Note that Xd = A 
when A < 1, since in that case b n = 1 and H = A. The 
asymptotic form for the distribution of the size of avalanches 
starting at node n, p n (x), can be obtained from the asymp- 
totic form of g n (s) as s — > 0. Therefore, we study Eq. (l32l by 
assuming g n (s) is small. In order to obtain an analytic expres- 
sion for the distribution of size we assume, in addition, that the 
network is close to critical, \\e> — 1| <C 1. Taking logarithms 
in Eq. d32l and using the approximation ln(l + g) ~ 3 — p 2 /2 
we obtain 

1 



9n(s) - ^g n {s) 2 



N 



N 



Hmn9m{s) H mn9m( s )- (33) 



As s and g n —> 0, the leading order terms are g n (s) = 
~ s + Y,m H mngm{s), or (H T - I)g = si, where g = 
[31,32, • ■ • ,g N ] T and 1 - [1, 1, ... , 1] T . When \X D -\\ « 1, 
and Ac is well separated from the rest of the spectrum of H, 
as we are assuming, g = s(H T — ~ v, where v is 

the left Perron-Frobenius eigenvalue of H (more precisely, 
we are assuming such a separation for A, but since H = A 
when A = 1 and we are assuming — 1| <C 1, by conti- 
nuity the assumption is valid for H as well). Since v is in- 
dependent of s, the solution up to first order is approximately 

9n{s) = g(s)v n /(v), where the term (v) = 5^=1"" is 
included to make g(s) independent of the normalization of v. 
For small s, and including the nonlinear terms, we expect the 
solution of Eq. ( 133b to be close to this solution, so we set 



9n{s) = g(s)(v) l v n 



(34) 



where e„ is a small unkown error term. Substituting Eq. 
into Eq. (l33l . using H T w = A_dv, and neglecting terms of 
order eg we get 

g&iv)-^ + e n {s) - \g{s) 2 {v)- 2 v 2 n 
= -s + A i g(s)(w) _1 i)„ 



N 



N 



H mn e m {s) ~ g{s) 2 {v)- 2 \ H ™v 2 m- (35) 



To eliminate the unknown error term e„, we multiply by the 
right eigenvector entry u n of H and sum over n. We use 
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FIG. 2: (color online) Histograms of avalanche duration shown above for networks of N = 10 5 nodes with power-law degree distribution, 
exponent 7 = 3.5 with Perron-Frobenius eigenvalues of A = 0.9 (left), A = 1.0 (center) and A = 1.1 (right). Symbols show the number 
of avalanches having duration d from a single simulation of 10 6 , 2 x 10 6 , and 10 6 avalanches, respectively, from left to right. Dashed lines 
provide a reference for the theoretical predictions described in Eqs. 10, {23), and l |13t . Note that the vertical position of the dashed lines 
was chosen arbitrarily. Due to predictions of exponential decay for the sub- and super-critical cases, the left and right plots are plotted on a 
log-linear scale, while the center plot is plotted on a log-log scale to show the power-law decay. Infinite duration avalanches in the supercritical 
case (right) are not displayed in the figure. 
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FIG. 3: (color online) Histograms of avalanche size shown above for networks of N = 10 nodes with power-law degree distribution, exponent 
7 = 3.5 with Perron-Frobenius eigenvalues of A = 0.9 (left), A = 1.0 (center) and A = 1.1 (right) on a log-log scale. Symbols show the 
number of avalanches having size x from a single simulation of 10 6 , 2 x 10 6 , and 10 6 avalanches, respectively, from left to right. Dashed 
lines provide a reference for the theoretical prediction x~ 3 ^ 2 exp(—x/x*) described in Eqs. i4l) and ((42). Note that the vertical position of 
the dashed lines was chosen arbitrarily. Infinite size avalanches in the supercritical case (right) are not represented in the data set. Agreement 
between theoretical prediction and measurement is excellent despite finite sample size noise. 



Hu = Adu and neglect (Ad — l)e„ to get 
= -s{u) + X D g(s)(v)~ 1 (uv) 



(36) 



where (xy) = j^^l n x n y n . Equation d36) is a quadratic 
equation for g(s), ag 2 + bg + c = 0, with 



2N(uv){v) 

b=(X D -l), 

(u)(v) 
c = -s- — -. 

{uv) 



(37) 
(38) 
(39) 



Solving for g(s) and substituting back into g n (s) = <\> n (s) — 1 
we find, choosing the root that guarantees g n < 0, 



1 



Asa- 



(uv) V n 



2a 



(v) 
(40) 



The moment generating function <j> n , first defined in Eq. 
can be interpreted as the Laplace transform of the distribution 
of size. Taking the inverse Laplace transform of the form of 
4> n (s) found in Eq. d40b we obtain that for large x, the distri- 
bution of size p n (x) is approximately given by 



p n (x) oc v n x 3 / 2 exp(-x/x*), 
where the characteristic size x* is given by 

(v)(u) 1 



4a- 



(41) 



(42) 
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The distribution of size is asymptotically an exponential times 
a power-law with exponent —3/2. Such a functional form de- 
scribes the distribution of the size of connected clusters near 
the percolation threshold in some network percolation models 
lfl9il20ll . In the critical case, when A = Xjj = 1, x* diverges 
and we recover a power-law distribution with exponent —3/2, 
which is the well-known exponent for critical branching pro- 
cesses ffl l29ll . It is interesting to note that this exponent, in 
our model, does not depend on the structure of the network, 
contrasting related percolation models where all nodes with 
the same degree are considered statistically equivalent 1I20I1 . 
Also note that the quantity a in Eq. (l42l depends implicitly on 
X D . 



V. NUMERICAL EXPERIMENTS 

In this section, we test the theoretical predictions of the pre- 
vious sections by directly simulating the process described in 
Sec.|n]on computer-generated networks. We first describe the 
processes used to construct networks and simulate avalanches. 

Networks were constructed in two steps. First, binary net- 
works (with adjacency matrix entries A mn G {0, 1}) were 
constructed via an implementation of the configuration model 
ll30ll . using N = 10 5 nodes, with nodal degrees drawn from 
a power-law distribution with exponent 3.5, i.e., the probabil- 
ity that a node has degree k is proportional to fc~ 3,5 . Second, 
each nonzero entry A mn was given a weight, drawn from a 
uniform distribution U[0, 1]. We then calculated the Perron- 
Frobenius eigenvalue of this weighted matrix, A, and mul- 
tiplied the matrix by A/ A, resulting in a matrix A with the 
desired eigenvalue A. We simulated avalanches for networks 
with A between 0.5 and 1.5, sampling more finely for values 
close to 1. 

Each simulated avalanche was created by first exciting a 
single network node, chosen uniformly at random, and then 
calculating the size and duration of the resulting avalanche as 
defined in Eqs. (T} and (ffji. If the resulting avalanche lasted 
for more than 10 6 time steps, we considered it as having in- 
finite duration and infinite size. In all cases, the initial ex- 
citation was included so that the minimum size was x = 1 
and the minimum duration was d = 1. For each subcritical 
(A < 1) and supercritical (A > 1) case, 10 6 avalanches were 
simulated, and for A = 1, we simulated 2 x 10 6 avalanches to 
better sample the very broad distribution of avalanche size at 
critic ality. 

A brief summary of the predictions of Sees. ITTTI and ITvlis 
as follows. The probability of an avalanche of duration d will 
decay as X d for subcritical networks (A < 1), as dr 2 for crit- 
ical networks (A = 1), and as A^, for supercritical networks 
(A > 1), where A^> is the Perron-Frobenius eigenvalue of the 
matrix D, Eq. (|3B- When |An - 1| < 1, the probability of 
a finite avalanche of size x will decay as a; -3 / 2 exp(— x/x*), 
where x* is a network-specific constant, given in Eq. (|42] >. 

In Figs.|2]and[3]we compare histograms of avalanche dura- 
tion and size obtained from direct numerical simulations for 
A = 0.9 (left), 1.0 (center), and 1.1 (right) with the theoret- 
ical predictions described in the previous paragraph (dashed 



lines). Note that, since our predictions allow for an un- 
specified proportionality constant, the vertical position of the 
dashed lines was chosen arbitrarily. In general, we find good 
agreement between the theoretical predictions of avalanche 
duration and size distributions with the histograms observed 
in the simulations. While the dashed lines in Figs. [2] and [3] are 
appealing to the eye, more quantitative measures of agreement 
between theory and experiment are shown in Figs. |4]and|5] 

To numerically test the agreement between theory and ex- 
periment for the distribution of avalanche duration, in Fig. [4] 
we compare the best fit A of the data to p(t) oc A*, calcu- 
lated through a nonlinear least-squares exponential regression 
on the simulated PDF of avalanche duration, to our theoretical 
predictions in Eqs. ((HJ and ( TBI (solid line). The agreement is 
excellent, though not exact, over the entire range of A values 
simulated. 

As a partial test of our theory for the distribution of 
avalanche size, we assume that the form of the distribution 
is x -3 / 2 exp(— x/x*), and estimate x* from the data, which 
we then compare with our theoretical prediction in Eq. (142V 
Noting that as A — > 1, x* will diverge, we estimated 1/x* 
via a nonlinear least-squares using Brent's minimization on 
the cumulative histogram of the avalanche size data. Since 
our theory describes only the asymptotic form of the distri- 
bution, this estimate was performed only on the largest 10% 
of measured data. [Similar results were obtained using the 
largest 5%, 1% and 0.1% of data (not shown), but when using 
more than the largest 10%, the minimizing x* value diverged, 
suggesting that we fit the power-law portion of data at the ex- 
pense of the exponential tail.] Figure shows the theoretical 
prediction (solid line) and the result of the numerical fit to the 
data (solid circles; the dashed lines are to aid the eye). As 
shown, agreement is quite good close to Xd = 1 (see the inset 



1.1 













Predicted 




• Measured 



0.4 0.6 0.8 1 1.2 1.4 1.6 

Perron-Frobenius Eigenvalue,^ 



FIG. 4: (Color online) A comparison of predicted duration decay 
rates [Eq. ^} and Eq. J13H (solid line), and numerical simulations 
(solid circles) plotted against A, the largest eigenvalue of the network 
adjacency matrix. Agreement is excellent for both the subcritical and 
supercritical numerical simulations. The distribution of avalanches 
durations decays as A* and A^ for A < 1 and A > 1, respectively, as 
indicated by arrows. 
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FIG. 5: (color online) Testing the prediction that avalanche size x is 
distributed as x~ 3 ^ 2 exp (— x/x*), we compare the theoretical pre- 
diction of x* (solid line) with x* estimated via regression on the 
largest 10% of avalanches from numerical simulations (solid circles, 
dashed line). Inset, identical data on a magnified domain around 
A = 1. Agreement is excellent for A near 1, and decreasingly accu- 
rate for much larger or smaller A. 



o 0.9- 



0.7 



0.6- 



■ Predicted 
Measured 



0.6 



0.8 1 1.2 

Perron-Frobenius eigenvalues 



1.4 



1.6 



FIG. 6: When the Perron-Frobenius eigenvalue A is larger than one, 
there is a non-zero probability of an avalanche starting at node n hav- 
ing infinite duration, as predicted by Eq. (|9}. Here we average the 
finite fraction of avalanches originating from node n over all nodes, 
showing excellent agreement between the fraction predicted by av- 
eraging Eq. ([9} (solid line) and fraction measured from simulation 
(solid circles). 



of Fig. 0, but less accurate for very subcritical or supercrit- 
ical networks. The latter is reasonable since the assumption 
that | Xd — 1 1 is a small quantity was used in the derivation of 
Eq. g2j. 

Although Figs. |4] and [5] demonstrate agreement between 
theory and measurement for supercritical networks, that anal- 
ysis was restricted to finite avalanches. To complement this re- 
sult, we compare the predicted fraction of infinite avalanches 
with the measured fraction, for various values of A^- The 
quantity b n in Eq. (O is the fraction of avalanches originat- 
ing at node n which will have finite duration and size. In 
Fig. [6] we show the fraction of avalanches that decay in fi- 
nite time, averaged over nodes, comparing theory (solid line) 
with experiment (solid circles). The theoretical fraction of 
avalanches was calculated by numerically solving Eq. (0 to 
find b n , n = 1 . . . , N, and then plotting Y^,n=i ^n/N as a 
function of A. The numerical fraction of finite avalanches was 
calculated by simulating 10 6 avalanches, each one starting at a 
random node (out of N = 10 5 nodes). If an avalanche lasted 
more than 10 6 steps, we counted it as an infinite avalanche. 
Then, an estimate of b n was calculated as the fraction of fi- 
nite avalanches starting at node n. The symbols in Fig. [6] 
show Yl n =i b n /N as a function of A. Agreement is excel- 
lent over the entire range of A values tested. Beyond aggre- 
gate statistics, we also test a more subtle prediction of Eq. (0. 
In Sec.|niJ we concluded that f n (t) = 1 — c n (t), the prob- 
ability that an avalanche started at node n lasts more than 
t steps, scales for large t as f n (t) oc A*w„, where v is the 
left Perron-Frobenius eigenvector of A. Other research in the 
network adjacency matrix literature has noted that the vector 
of nodal out-degrees (in-degrees) is a good approximation for 
the right (left) dominant eigenvector of A in the absence of 
degree correlations ||3lll . In this light, our prediction above 



is understandable: when there are not degree correlations in 
the network, a node with a larger right eigenvector entry (and 
thus larger out-degree) will tend to produce longer avalanches. 
Therefore, in order to fully test our prediction, we created net- 
works with assortative mixing by degree l32ll . a type of degree 
correlation which we measure using the coefficient p ||3lll . 

I uinh,out\ 

" — /Uin\ lhout\ ' ^ ' 

\ K n hX^m h 

where (-) e denotes an average over all edges and k are 
weighted nodal degrees defined as fc" 1 = J2 m A m n an d 
kn Ut = Em A nm . In the absence of degree correlations 
between connected nodes (fc^fc^ 1 *),, = (k™) e (^m*) e an d 
p = 1. In assortative networks, there exists a positive cor- 
relation (p > 1) between the in-degree at node n and the 
out-degree at node m at the ends of a directed link from n 
to m. When the correlation is negative {p < 1), the network 
is called disassortative. Thus we created Erdos-Renyi random 
networks with N = 10 4 nodes, and rewired each network via 
a link-swapping process (as described in Ref. 1I3TI0 until we 
had very assortative and disassortative networks (p = 1.2 and 
p = 0.8, respectively). Eq. (|7]i implies that in such networks, 
the tails of the CDF of avalanches originating at node n will be 
proportional to the corresponding entry of the right eigenvec- 
tor, which may differ significantly from the nodal out-degree. 
For the a subcritical network (A = 0.95) with assortativity 
coefficient p = 0.8 we plot / n (30) and its corresponding en- 
try in the right dominant eigenvector v n for each node n, in 
Fig. [7] showing that proportionality is excellent. In the inset 
of the same figure we plot /„(30) against the corresponding 
out-degree k° ut for each node n, showing that proportionality 
to out-degree does not hold. Assortative networks produce the 
same effect, but are not shown here. 
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FIG. 



right dominant eigenvector entry v n 

7: Testing the node-specific prediction of Eq. (0, avalanches 



were simulated on a subcritical (A = 0.95) and disassortative (p = 
0.8) Erdos Renyi random network with TV = 10 4 nodes. In the large 
plot, the fraction of avalanches originating at node n that last longer 
than 30 time steps, /„(30), is plotted against the corresponding en- 
try in the right Perron-Frobenius eigenvector, v n . In the inset, the 
same values /„(30) are plotted against the corresponding out-degree 
The eigenvector entry v n does a significantly better job than 
out-degree fc°j Ui of predicting the duration of avalanches originating 
at node n in disassortative networks (shown) and for assortative net- 
works (not shown). 



VI. DISCUSSION 

We have presented an analysis of the asymptotic distribu- 
tions of the duration and size of avalanches in complex net- 
works. This work is of interest in various applications, most 
notably neuroscience lRl- [l3ll and the analysis of power-grid 
failure cascades JUj]. While some of our results, such as the 
functional forms for the distributions, are analogous to those 
found in classical Galton- Watson branching processes lfl4h or 
in mean-field models lEoll . we emphasize the distinguishing 
aspects of our results: (i) We generalize the criterion for criti- 
cality to A = 1, which depends on the topology of the network 
in ways that previous results do not capture. For example, in 
critical branching processes ll29l the condition for critic ality 
is (d) = 1. (ii) The parameters of the asymptotic distributions 
in the various regimes are affected by the network topology, 
and our results allow us to predict how various factors such as 
network degree distributions or degree-degree correlations af- 
fect these parameters [e.g., the parameter x* in Eq. (l42l or Ad 
in Eq. (fTJtl. (iii) In contrast to previous studies, our results 
allow us to predict the statistics of avalanches generated at a 
particular node. This might be of critical importance in cer- 
tain applications where the adjacency matrix is known or can 
be inferred (such as the power grid or the Autonomous System 
network of the internet) since one can then allocate resources 
to prevent avalanches, if so desired, that start at the nodes 
which tend to generate the largest avalanches. As shown in 
Fig. |7] the naive prediction that the nodes with the largest out- 
degree generate the largest avalanches is not necessarily true 
when the networks have nontrivial structure, such as degree 



correlations. 

In developing our theory, we made some assumptions 
which we now discuss. First, we assumed that the network 
was locally tree-like. This allowed us to treat avalanches prop- 
agating to the neighbors of a given node as independent of 
each other. While this is a good approximation for the net- 
works we used, it is certainly not true in general. In particular, 
avalanches propagating separately from a given node might 
excite the same node as they grow. The result is that the num- 
ber of nodes that the avalanches excite in the simulation may 
be less than what the theory would predict. In running our 
simulations, we addressed this issue in two ways: first, we 
kept track of the number of times two branches of the same 
avalanche simultaneously excited the same node n, finding 
it to be an increasing function of avalanche size and Perron- 
Frobenius eigenvalue, yet still negligible when compared to 
the total number of excitations. In addition, each time such an 
event occurred, we separately generated an avalanche start- 
ing from the doubly excited node n and corrected both the 
size and duration of the original avalanche by incorporating 
these additional avalanches. We found that doing this had 
no appreciable effect on the measured distributions, and so 
all figures shown in this manuscript are produced from sim- 
ulation data without the additional compensating avalanches 
included. This, and the fact that the numerical simulations 
are described well by the theory, suggest that the interaction 
of avalanches propagating to different neighbor nodes can be 
safely neglected in the networks studied. The performance of 
our theory in networks that are not locally tree-like, such as 
networks with a high degree of clustering, is left for future 
research. Another approximation we used is that the Perron- 
Frobenius eigenvalue A is well separated from the rest of the 
spectrum. This is a good approximation in networks without 
well defined communities, but can break down in networks 
with strong community structure lf27ll . 

Finally, we note that our results show that the experimental 
signatures of criticality in neural systems (characterized by 
a power-law distribution of avalanche size and duration with 
exponents —3/2 and —2, respectively fllUliilliltl) are robust 
to complex underlying network topologies. 
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Shew. DBL and MC were supported by NSF MCTP Grant No 
DMS-0602284. JGR was supported by NSF Grant No. DMS- 
090822. EO was supported by ONR MURI Gran N00014-07- 
1-0734. 



Appendix A: Probability of finite avalanche duration 

In this Appendix we establish that the probability of finite 
avalanches, under our assumptions, is always one when A < 1 
(critical and subcritical networks), and becomes less than one 
when A > 1 (supercritical networks). These probabilities, 
b n = linif^oo c n (t), satisfy the equation 



N 

n[< 



1 A mn ) + A Tl 



(Al) 
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First, we show that if A < 1, where A is the Perron-Frobenius 
eigenvalue of A, then the only solution to the equation above 
is b n = 1. Letting b n = !-/„, we have for all n 



N 

n[ 



1 Amn fn 



l -In 



Using the Weierstrass product inequality 113 311 . 



(A2) 



N 

E 

m — 1 



Amnfm — 1 



N 

n[ 



1 A mn fjji — fn 5 



(A3) 



with equality only if (i) A mn f m = for all m, or (ii) 
Ajm/m = for all m ^ k and A kn f k = 1 for some k 
113311 - If u is the right Perron-Frobenius eigenvector of A, this 
implies, since u T A T = Au T , 



i T A T { = Au T f > u T f. 



(A4) 



If there is a nonzero f n , then u T f > since the Perron- 
Frobenius eigenvector has positive entries for irreducible A. 
Therefore, if A < 1 we must have /„ = for all n. If A = 1, 
Eq. dA4l ) implies equality in Eq. dA3b . which implies either 
(i) A mn f m = for all m, and thus /„ = by dA31 >. or (ii) 
A mn fm = for all m ^ k and A kn fk = 1 for some k, 
which is impossible since we assumed that the entries of A 
are strictly less than one and f k is a probability. Therefore, 
we must have /„=0ifA = l,a valid argument for any n. 
Together with the previous argument above, we conclude that 
b n = 1 for all 77. if A < 1. 



Now, we show that if A > 1 then lim t _ i . 00 c n (t) = b n < 1. 
To show this, we view Eq. (HJi as a dynamical system, and note 
that the analysis of Sec. lIII Al applied to the case A > 1, shows 
that the fixed point b n = 1 is linearly unstable. If we show that 
Cn(t) is nondecreasing with t, then the limit b n must be less 
than one. By induction, we will prove that c„(i + 1) > c„(i) 
for all n. First, we have c„(0) = and c„(l) = [] m (l — 
Anm) > 0, so the statement is valid for t = 0. Then, assume 
c m {t) > c m {t — 1) for all m and consider c n {t + l)/c n (t), 
noting that c n (t) > 0: 



JY 



C n (t + 1) _ T7 (1 Amn) + Amn c m(t) 
Cn(t) 



n 5 



n 



j (1 Amn) + A m n c m{t 1) 

A m n{Cm(t) ~ C m {t ~ 1)) 



N r 



m— 1 
> 1, 



(l A mn j -\- A mn c m (t l) 



(A5) 



which proves the desired statement. Note that, although from 
the definition (O, it follows that c n (t) are nondecreasing, this 
proof is necessary since Eq. <j4j is an approximation. 



Appendix B: A > 1 => \d < 1 

In this Appendix we argue that the Perron-Frobenius eigen- 
value of the similar matrices H and D is less than one when 
the Perron-Frobenius eigenvalue of A is greater than one: 
A > 1 =>■ Ad < 1. Recall that the matrix D was defined 

as 



Dmn. — 



h A 

Vn-f^-mn 



(1 A ra n) ~\~ bmA Tl 



(Bl) 



where b n , the probability that an avalanche starting at node n 
is finite, satisfies 



b n = ] — Amn) + A m nbr> 



(B2) 



Now, suppose that A is such that A > 1, and introduce a 
parameter a < 1 by defining b n (a) as the b n corresponding 
to the matrix aA, which satisfies 



b n {a) = ] [(1 -aAmn) +o>A mn b m (a) 



(B3) 



Now, calculate the derivative of b n (a) with respect to a, 



dbnja) . / n 
- b n {a) 

m=l 
db n 



-A mn + Amnbm(a) + aA 



db m (a) 
™ da 



da 



(1 - aAmn) + aAmnb m {a) 



(B4) 

Letting \i n = =1 > and evaluating the expression above 
at a = 1, we get 



-A r , 



Amnf^i] 



( 1 A m .n ) ~1~ A m n &n 



N 
m—1 

N N 

- ^ D rnn (pm 1) ~\~ ^ D mn fJ, rn . 
m—1 m—1 



In matrix form, 



(D T - I)n = D T (1 - b), 



(B5) 
(B6) 

(B7) 



where 1 = [1,1,..., 1] T , b = [b%, 62, • • • , &at] t , and \i = 
[fj,i, fi2, ■ ■ ■ j Mjv] T - Now, we left multiply the previous equa- 
tion by u T , where u is the right Perron-Frobenius eigenvector 
of D, satisfying u T D T = AdU T , to get 



(X D - l)u T n = Au T (l -b) 



(B8) 



If A > 1, AppendixlAlshows that the entries of (1 — b) are all 
positive. Since the Perron-Frobenius eigenvector u has posi- 
tive entries as well (since we are assuming A is irreducible), 
the right hand side of Eq. dB8b is positive. Now, we argue 
that the vector /i has nonpositive entries: as a increases, the 
probability of an excitation passing between any pair of nodes 
increases, and thus the probability of having a finite avalanche 
can not increase, i.e., db n /da < 0. Therefore, the term u T /i 
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on the left hand side must be nonpositive and, since the right — 1 must be negative, that is, < 1, 

hand side is nonzero, it must be negative. Thus, the term 
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